Getting started

Adapted from: SOSTOOLS' SOSDEMO1 (See Section 4.1 of SOSTOOLS User's Manual) and Example 2.4 of [PJ08]

P. Parrilo and A. Jadbabaie Approximation of the joint spectral radius using sum of squares. Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402

using DynamicPolynomials
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

\[ 5y^{4} - x^{2}y^{2} + 2x^{3}y + 2x^{4} \]

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
con_ref = @constraint(model, p >= 0)
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We see above that the solver found a feasible solution. We now inspect this solution:

q = gram_matrix(con_ref)

GramMatrix with row/column basis:
 SubBasis{Monomial}([y^2, x*y, x^2])
And entries in a 3×3 SymMatrix{Float64}:
  5.0                     4.4408920985006264e-17  -2.9518518518518517
  4.4408920985006264e-17  4.903703703703703        1.0
 -2.9518518518518517      1.0                      2.0

We can get the SOS decomposition from the gram matrix as follows:

sosdec = SOSDecomposition(q)

(-2.1236350925148186*y^2 + 0.685602733872361*x*y + 1.404328747693391*x^2)^2 + (-0.6931517747987985*y^2 - 2.1053486943394364*x*y - 0.020343251432301153*x^2)^2 + (0.0985627258797928*y^2 - 0.034050994900062616*x*y + 0.16567112157245356*x^2)^2

We now seek for the SOS decomposition of the following polynomial:

p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

\[ y^{2} + x^{2} - xy^{2} + 4x^{4}y^{6} \]

We build the same model as previously with this new polynomial. Here we can use Model instead of SOSModel as we explicitly constrain p to belong to the SOS cone with p in SOSCone().

model = Model(solver)
con_ref = @constraint(model, p in SOSCone())
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can query the SOS decomposition directly from the constraint reference as follows:

sos_decomposition(con_ref)

(-0.8633544111907671*y + 0.19297758710767499*x*y + 3.0514841249382645e-18*x*y^2 + 1.9748040632461146*x^2*y^3)^2 + (-1.8744116508952374e-16*y + 0.7962855107812752*x - 1.0528501085696115e-16*x*y - 1.7452642425911875*x*y^2 - 8.635065971616514e-17*x^2*y^3)^2 + (0.24383463418983628*y - 1.7639164207923265e-16*x - 1.548639212897575*x*y + 3.741600489401518e-17*x*y^2 + 0.25793362243589013*x^2*y^3)^2 + (-5.3930412033044254e-17*y + 0.6049209744419549*x - 7.96731030959412e-17*x*y + 0.27599821010522496*x*y^2 - 9.028384658229797e-18*x^2*y^3)^2 + (-0.4417735074073057*y - 1.0850742441349411e-16*x - 0.10009637589813172*x*y + 1.468808219015887e-17*x*y^2 - 0.18335527863613862*x^2*y^3)^2

This page was generated using Literate.jl.